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ABSTRACT 


Previous  methods  of  determining  the  wellbore 
temperature  distribution  have  involved  a  compromise  between 
accuracy,  and  the  time  required  to  obtain  the  information. 

To  overcome  this  problem,  a  computer  model  has  been 
developed  which  utilizes  a  direct  solution  technique  to 
solve  the  finite  difference  equations  describing  the 
transient  heat  transfer  mechanisms  in  a  wellbore  during 
drilling  operations.  The  method  involves  representing  this 
system  of  equations  by  a  heptadi agona 1 ,  asymmetric  band 
matrix,  and  solving  this  matrix  by  means  of  a  band 
algor i thm. 

A  significant  improvement  in  the  solution  time  is 
achieved,  thus  allowing  the  model  to  be  used  at  the  wellsite 
to  provide  dynamic  wellbore  temperature  distributions. 

The  results  of  a  parametric  sensitivity  analysis 
carried  out  using  the  model  indicate  that  a  number  of 
assumptions  made  in  previous  models  are  invalid,  and  that 
certain  factors  ignored  in  previous  models,  have  a 
significant  effect  on  the  wellbore  temperature  distribution. 
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1.  INTRODUCTION 


Today  the  drilling  of  a  well  to  depths  approaching 
30,000  feet  on  land,  and  offshore  in  water  depths  of  several 
thousands  of  feet,  is  becoming  increasingly  commonplace. 
This,  coupled  with  the  fact  that  the  search  for  oil  is  being 
concentrated  in  the  hostile  environments  of  such  areas  as 
the  North  Sea  and  the  High  Arctic,  has  made  the  drilling  of 
a  well  a  highly  complicated  and  expensive  operation. 

With  such  enormous  costs  and  the  highly  advanced 
technology  involved,  it  is  necessary  that  the  most  efficient 
means  possible  be  employed  in  the  drilling  of  oil  wells. 
Hence,  accurate  procedures  for  estimating  the  numerous 
variables  involved  are  becoming  essential,  both  at  the 
design  stage  and  at  the  wellsite,  even  though  the 
complicated  nature  of  some  of  these  variables  would  mean  the 
development  of  intricate  models. 

Although  various  models  have  been  formulated  to  enable 
certain  parameters  to  be  estimated,  most  have  been 
oversimplified  so  as  to  enable  their  use  at  the  wellsite  in 
the  form  of  simple  equations,  solvable  by  a  sliderule, 
nomographs  or  plots.  As  a  result  the  practical  application 
of  the  relevant  theory  is  minimal.  In  the  present  era  of 
powerful  minicomputer  systems,  this  can  no  longer  be 
justified  and  it  is  becoming  increasingly  more  feasible  to 
utilize  the  complicated  theory  that  exists  in  the 
literature,  at  the  wellsite.  Models  developed  from  such 
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theory,  with  subsequent  improvements  in  accuracy,  are 
probably  most  useful  at  the  well  planning  and  design  stage 
in  view  of  the  fact  that  large  powerful  computers  are 
available  with  virtually  no  limitations  on  solution  time  or 
storage  space,  for  all  practical  purposes. 

One  of  the  parameters,  of  which  a  thorough  Knowledge 
would  be  very  helpful,  is  temperature.  In  the  past  it  has 
been  convenient  to  ignore  the  temperature  distribution  in  a 
well  and  to  assume  an  isothermal  system,  primarily  because 
no  practical  means  existed  for  determining  the  wellbore 
temperature  profile.  An  accurate  means  of  estimating  the 
temperature  distribution,  and  its  variation  with  time,  would 
have  a  variety  of  applications  as  follows: 

1.  Enable  the  dynamic  temperature  profile  and  bottomhole 
temperature  to  be  determined  rather  than  the  static 
temperatures  presently  available  from  electric  logging 
tools . 

2.  Improve  cementing  programme  design,  particularly  with 
regard  to  the  amount  of  retarder  required,  and  the 
setting  time. 

3.  Improve  drilling  fluid  design  by  providing  information 
on  the  actual  circulating  temperatures  so  as  to  enable 
high  temperature  modifications  to  be  made  to  the 
drilling  fluid  programme. 

4.  Enable  casing  thermal  stresses  to  be  determined. 

5.  Provide  improved  well  design  in  permafrost  regions. 
Improve  injection  and  production  operations. 


6. 


3 


At  the  present  time  the  only  temperature  information 
available  at  the  wellsite  is  the  inlet  and  outlet  fluid 
temperatures  and  the  static  bottomhole  temperature  obtained 
from  a  logging  tool.  The  purpose  of  this  dissertation  is  to 
employ  the  existing  theory  in  the  literature  on  transient 
heat  flow  in  a  wellbore  to  develop  a  computer  model  to 
estimate  the  temperature  distribution  and  its  variation  with 
circulating  time.  Unlike  previous  models,  it  will  be 
feasible  to  use  this  model  at  a  remote  wellsite  without  a 
significant  loss  of  accuracy.  Furthermore  this  model 
incorporates  a  number  of  improvements  over  previous  models 
in  terms  of  both  accuracy  and  applicability. 

A  further  aim  of  this  study  was  to  use  the  computer 
model  to  simulate  wellbore  temperature  distributions  so  that 
a  complete  parametric  sensitivity  analysis  could  be  carried 
out  to  indicate  which  parameters  have  the  most  significant 
effect  on  the  temperature. 
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2.  REVIEW  OF  THE  AVAILABLE  MODELS 


Temperature  distribution  in  a  wellbore  is  a 
complicated  function  of  a  large  number  of  variables,  many  of 
which  are  not  Known  with  any  degree  of  certainty.  For  this 
reason,  any  mathematical  model  must  contain  simplifying 
assumptions  to  enable  a  solution  to  be  obtained.  The  more 
complicated  the  model,  the  more  data  are  required.  In  view 
of  the  fact  that  any  model  can  only  be  as  accurate  as  the 
data  that  it  uses,  the  intricacy  of  the  model  is  not  in 
itself  a  sufficient  criterion  for  accuracy. 

A  number  of  models  exist  for  estimating  bottomhole 
temperatures  and  wellbore  temperature  distributions  in  the 
literature.  Some  of  these  have  been  designed  for  practical 
application,  and  each  includes  a  variety  of  simplifying 
assumptions.  The  earlier,  hand-solved  methods  tend  to  be  the 
most  simplified  and  inaccurate. 

The  bottomhole  temperature  charts  developed  by  Farris 
(1941)  were  based  solely  on  measurements  of  the  circulating 
temperatures  in  five  shallow  Gulf  Coast  wells.  Despite  the 
obvious  inaccuracies  and  oversimplifications  involved  in  the 
use  of  the  Farris  charts,  and  their  restrictive  application, 
the  API  recommends  them  for  determining  setting  schedules 
for  cement  slurries.  The  severe  shortcomings  of  these  charts 
prompted  research  into  a  more  accurate,  mathematical  method 
for  estimating  circulating  temperatures . 
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Edwardson  et  al.  (1962)  solved  the  differential  heat 
conduction  equations  to  estimate  the  formation  temperature 
distribution  before  and  after  circulation.  However,  this 
work  was  more  concerned  with  temperature  distribution  in  the 
formation  as  a  result  of  drilling  fluid  circulation,  and 
does  not  allow  the  direct  calculation  of  the  wellbore 
temperature  distribution.  Moreover,  assumed  formation 
temperature  profiles  were  used  at  the  end  of  circulation. 

Tragesser  et  al.  (1969)  expanded  Edwardson' s  method  so 
as  to  allow  for  such  variables  as  depth,  pump  rate,  hole 
size  and  so  forth,  but  their  method  involved  the  same 
limiting  assumptions. 

Holmes  and  Swift  (1970)  assumed  steady  state  conditions 
in  the  wellbore  and  surrounding  formation,  and  solved  the 
relevant  heat  transfer  equations  analytically,  using  this 
assumption.  This  method  has  the  advantage  of  being  simple 
and  more  accurate  than  previous  methods,  but  the  assumption 
of  steady  state  heat  flow  is  a  critical  one  which  is  only 
satisfied  after  impract ical ly  long  circulation  times. 

In  terms  of  the  theoretical  development  of  the  heat 
transfer  relationships  for  a  wellbore,  the  work  of  Raymond 
(1969)  is  the  most  comprehensive.  It  is  this  work  that  has 
served  as  a  basis  for  all  the  recent  research  on  the 
subject.  He  advocated  the  use  of  the  principle  of 
superposition  and  the  Hurst  and  van  Everdingen  functions  to 
solve  numerically  for  unsteady  state  conditions.  However,  he 
contended  that  the  pseudo-steady  state  condition  was 
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accurate  enough  for  all  practical  purposes  and  described  an 
analytical  means  of  solving  for  this  condition. 

Most  recent  research  has  dealt  with  means  of  solving 
the  unsteady  state  equations  formulated  by  Raymond,  or 
variations  on  these,  using  finite  difference  techniques, 
e.g.  Sullivan  (1970),  Keller  et  al .  (1973),  Sump  and 
Williams  (1973)  and  Oster  (1976).  The  approach  of  Keller  et 
al.  is  the  most  comprehensive  of  these  methods.  As  well  as 
solving  the  finite  difference  equations,  variations  in  the 
wellbore  geometry  were  included.  A  further  refinement  was 
the  allowance  for  heat  generation  within  the  fluid  column 
due  to  friction  forces  and  the  rotational  energy. 

It  is  the  work  of  Keller  et  al.,  which  itself  is  an 
extension  of  Raymond's  work,  that  is  the  basis  of  this 
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3.  STATEMENT  OF  THE  PROBLEM 


A  review  of  the  literature  on  determining  static  and 
dynamic  wellbore  temperature  distributions  suggests  that  a 
compromise  must  be  made  between  accuracy  and  time  required 
to  obtain  the  temperature  distribution.  The  early  methods 
are  inaccurate  but  the  temperature  data  is  readily  acquired; 
the  more  recent  ones  are  quite  accurate  but  necessitate 
impract ical ly  long  periods  of  time  to  obtain  solutions. 

A  need  exists  for  a  method  that  is  both  accurate  and 
fast  thus  enabling  it  to  be  used  at  the  wellsite  where  it 
would  be  of  greatest  value.  Such  an  application  is  not 
possible  with  any  of  the  existing  means  of  obtaining  dynamic 
temperature  data,  without  sacrificing  accuracy. 

This  dissertation  utilizes  the  existing  theory  on  the 
subject  of  wellbore  transient  heat  transfer  to  develop  a 
computer  model  that  is  both  accurate,  and  by  means  of  a 
direct  solution  technique,  much  faster  than  previous 
mathematical  methods. 

The  study  also  investigated  the  effect  of  considering 
the  non-Newtonian  flow  behaviour  of  the  circulating  medium 
on  the  temperature  profiles  generated  by  the  model. 
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4.  THEORY  AND  DEVELOPMENT  OF  THE  MODEL 


The  physical  system,  upon  which  this  computer  model  is 
based,  is  shown  schematically  in  Fig.  4-1 

The  flow  of  a  drilling  fluid  in  a  wellbore  may  be 
divided  into  three  distinct  regions  as  follows: 

1.  Downward  flow  through  the  drill  string. 

2.  Flow  from  the  drill  string  through  the  drill  bit  into 
the  annulus. 

3.  Upward  flow  through  the  annulus. 

The  temperature  of  the  fluid  in  each  region  is  dependent 
upon  a  number  of  different  thermal  processes.  The  fluid 
enters  the  drill  pipe  with  a  Known  temperature  and  its 
change  in  temperature  is  determined  by  the  rate  of  thermal 
convection  down  the  fluid  column  and  the  rate  of  convective 
heat  transfer  radially  between  the  fluid,  the  pipe  wall  and 
the  annulus.  Vertical  and  radial  heat  conduction  within  the 
pipe  wall  is  also  present. 

As  the  fluid  flows  up  the  annulus,  its  temperature  is 
dependent  upon  the  rate  of  heat  convection  up  the  annulus, 
the  rate  of  radial  convection  between  the  annulus  fluid,  the 
drill  pipe  wall  and  the  fluid  within  the  drill  pipe,  the 
rate  of  radial  convection  between  the  annulus  fluid  and  the 
formation  or  casing  and  the  conduction  through  the 
formation . 

Since  transient  heat  transfer  is  considered,  the  time 
of  circulation  has  an  important  effect  on  the  temperature  of 
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FIG. 4-1 : 


SCHEMATIC  DIAGRAM  OF  WELL- 
FOR  DIMENSIONS  SEE  TABLE  4-1 
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the  fluid  as  has  the  heat  generation  within  the  three  flow 
regions  due  to  frictional  forces  and  the  rotational  energy 
of  the  drill  string  and  the  drill  bit. 


4.1  Assumptions 

To  arrive  at  the  energy  balances  describing  the  thermal 
behaviour  of  the  wellbore  certain  assumptions  about  the  heat 
transfer  mechanisms  and  flow  behaviour  are  required. 

1.  Flow  is  steady  state  and  fully  developed. 

2.  Flow  is  turbulent  in  the  drill  pipe  and  drill  bit  and 
laminar  in  the  annulus. 

3.  Heat  transfer  within  the  drilling  fluid  is  by  axial 
convection.  Axial  conduction  may  be  neglected. 

4.  The  radial  temperature  gradient  within  the  drilling 
fluid  may  be  neglected. 

5.  Heat  generation  by  viscous  dissipation  within  the  fluid 
may  be  neglected. 

6.  Fluid  properties  such  as  density,  thermal  conductivity 
and  heat  capacity,  are  independent  of  temperature. 

4.2  Energy  Balances 

The  partial  differential  equations  describing  the 
energy  balances  within  the  system  are  derived  in  Appendix  I. 
Equation  4.1  describes  the  heat  flow  within  the  drill 
string. 
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The  terms  on  the  left  represent  the  vertical  and  radial 
convective  heat  transfer,  respectively.  The  right  hand  terms 
represent  the  accumulation  of  energy  within  the  drill  string 
and  the  energy  sources,  respectively. 

The  energy  balance  in  the  drill  pipe  wall  is  given  by 
Equation  4.2. 


a2T 


w 


2r  h 
P  P 


2r  h 
a  a 


w 


3Z‘ 


-  (T  -  T  )  +  - 2_2_  (t 

(r2  -  r2)  p  w  (r2  -  r2)  a 


3T 


=  p  C 

p  w 


w 


at 


(4.2) 


The  first  term  on  the  left  hand  side  represents  the  vertical 
heat  conduction  within  the  wall.  The  second  and  third  terms 
account  for  the  radial  heat  transfer  to  the  drilling  fluid 
inside  and  outside  the  wall.  The  right  hand  term  represents 
the  heat  accumulation. 

The  energy  balance  in  the  flow  annulus  is  expressed  by 
Equation  4.3. 
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The  three  left-hand  terms  represent  the  vertical  convective 


heat  transfer  within  the  fluid,  radial  convection  between 


the  drilling  fluid  and  the  drill  pipe  wall  and  radial 


convection  between  the  drilling  fluid  and  the  casing  or 
formation,  respectively.  Heat  accumulation  and  generation 
are  accounted  for  by  the  two  right-hand  terms. 

Equation  4.4  is  a  two-dimensional  thermal  conductivity 


f  =  pfcpf  3Tf 


(4.4) 


equation  representing  heat  flow  in  the  formation.  The  two 
left-hand  terms  account  for  the  vertical  and  radial 
conduction,  and  the  right-hand  term  accounts  for  the  heat 
accumulation . 

Since  the  original,  theoretical  work  of  Raymond  (1969) 
on  this  subject,  these  equations  have  formed  the  basis  of 
the  various  thermal  wellbore  models  that  have  been 
developed.  Prior  to  Raymond's  work  the  thermal  models  were 
all  very  simplified  to  the  extent  that  they  were  of  little 
value  in  estimating  the  temperature  distribution  with  any 
degree  of  accuracy. 

4.3  Solution  of  the  Partial  Differential  Equations. 

The  solution  of  these  equations  to  obtain  a  temperature 
distribution  as  a  function  of  time  is  complicated  and  the 
models  developed  subsequent  to  Raymond's  paper  incorporated 
solution  methods  based  on  finite  difference  techniques.  The 
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solution  of  these  finite  difference  equations  all  involved 
iterative  methods  with  the  obvious  disadvantages  of  long 
solution  times  and  the  concurrent  problems  of  stability  and 
accuracy.  Keller  et  al .  (1973),  whose  work  on  the  subject  is 
the  most  comprehensive  to  date,  quote  a  solution  time  of  170 
seconds  on  a  CDC  6400  for  a  sample  problem.  Obviously  such  a 
method  is  impractical  for  any  wellsite  application. 

In  the  present  solution  method,  the  wellbore  and  the 
adjacent  formation  are  represented  by  a  two-dimensional, 
mesh-centred  grid  consisting  of  ten  radial  elements  and  a 
variable  number  of  vertical  elements  depending  on  the  well 
depth,  but  approximately  one  element  for  each  100  feet  of 
depth.  Each  of  the  radial  elements  corresponds  to  a 
different  portion  of  the  wellbore  cross-section  from  the 
centre  of  the  drill  string  into  the  formation,  as  shown  in 
Fig.  4-2.  The  ten  radial  elements  allow  for  a  well  profile 
with  up  to  three  casing  strings  and  the  resulting 
mesh-centred  grid  is  of  a  form  similar  to  that  in  Fig.  4-2. 
Using  this  grid  as  a  basis,  the  above  partial  differential 
equations  can  be  written  in  finite  difference  form  for  each 
element  of  the  grid  so  as  to  describe  the  transient  heat 
flow  in  that  element.  The  finite  difference  equations  are 
written  in  an  implicit  form  and  are  derived  in  their  general 
form  in  Appendix  II.  Two-point  forward  and  two-point 
backward  difference  approximations  are  used  to  represent  the 
first  order  spatial  derivatives  and  the  time  derivatives. 

The  second  order  spatial  derivatives  are  represented  by 
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FIG. 4-2:  MESH  CENTRED  GRID  REPRESENTING 
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three-point  central  difference  approximations. 

The  finite  difference  equations  may  be  arranged  into 
the  following  general  form  for  each  node. 


n+1  n+1  n+1  n+1  n+1 

B..T.  .  i  +  D..T.  ,  .  +  E..T.  .  +  F.  .T. ,,  .  +  H..T. 
ij  i,J-l  ij  i-l »J  iJ  i,J  i+l  ,J  iJ  i,J+l 


(4.5) 


Using  a  standard  ordering  procedure  the  coefficients  in  the 
above  general  equation  for  each  node  are  represented  by  a 
pentadi agona 1  matrix,  so  that  the  original  N  x  10  block 
matrix  becomes  a  ION  x  ION  band  matrix  of  the  form  shown  in 
Fig.  4-3. 

Having  set  up  the  band  matrix,  certain  modifications  to 
it  are  required  to  satisfy  the  initial  and  boundary 
conditions.  These  conditions  are  shown  schematically  in  Fig. 
4-4. 

Equation  4.6  describes  the  initial  temperature 

T1.  .  =  Ts  +  mz,  (4.6) 

i  ,J  J 

profile  as  being  the  geothermal  gradient. 

The  boundary  condition  at  the  upper  surface  is  given  by 


aT 

az 


z=0 


0 


(4.7) 


Equation  4.7,  which  indicates  no  flow  of  heat  between  the 
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FROM  3  X  10  BLOCK  MATRIX  USING 
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FIG. 4-4:  INITIAL  AND  BOUNDARY  CONDITIONS 
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earth  and  the  atmosphere. 

The  inner  and  outer  radial  boundary  conditions  are 
described  by  Equations  4.8  and  4.9. 
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These  equations  imply  that  no  heat  flows  across  these 
boundar i es . 

Equation  4.10  denotes  the  inlet  fluid  temperature  as 
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n 

1,0 


(4.10) 


being  the  boundary  condition  at  the  innermost  upper  element. 

The  boundary  condition  along  the  lower  surface  is  given 
i n  Equation  4.11. 
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(4.11) 


Equation  4.12  ensures  that  the  temperatures  are 
equalized  at  the  drill  bit. 


n 

T,  • 

1  ,jmx 


n  n 

T2,jmx  ~  T3,jmx 


(4.12) 


Because  of  this  final  condition,  the  band  matrix  must  be 
modified  slightly  by  "folding  in"  the  H  coefficients  of  the 
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second  and  third  radial  elements  of  the  final  row.  To 
achieve  this,  two  extra  lower  co-diagonals  are  required, 
each  with  one  of  the  coefficients.  The  matrix  in  its  final 
form,  incorporating  the  initial  and  final  boundary 
conditions,  is  an  asymmetric,  heptadi agonal  band  matrix. 

Because  of  the  potentially  large  storage  requirements 
of  such  a  matrix,  especially  where  a  deep  well  is  being 
modeled,  it  is  transformed  into  a  band  storage  matrix  of  the 
form  shown  in  Fig.  4-5,  thus  reducing  its  size  from  ION  x 
ION  to  ION  x  23.  This  band  storage  matrix  is  then  solved  by 
a  direct  band  algorithm  technique  involving  factorization 
into  lower  and  upper  triangular  matrices.  The  details  of 
this  method  are  adequately  covered  in  the  literature  e.g. 
Peaceman  (1973),  and  will  not  be  considered  here.  The  actual 
solution  was  achieved  by  means  of  an  IMSL  subroutine  LEQT1B, 
and  is  described  in  more  detail  in  Appendix  III. 

When  the  band  algorithm,  direct  solution  method  was 
used  for  the  test  data  in  Table  4-1,  the  resulting  annulus 
temperature  profile  was  identical  to  that  of  Keller  et  al . 
(1973).  However,  instead  of  requiring  a  solution  time  of  170 
seconds  on  a  CDC  6400,  about  3  seconds  were  needed  on  an 
Amdahl  470. 

4.4  The  Selection  of  an  Appropriate  Rheological  Model 

The  basis  of  any  model  of  the  drilling  process  is  the 
drilling  fluid  or  mud.  The  drilling  fluid  serves  a  variety 
of  purposes  such  as  cooling  the  bit,  transporting  formation 


'  •• 


■  ’  ■'  :  :  ;  l 


' 


20 


FIG. 4-5:  30  X  21  BAND  STORAGE  MATRIX 

FROM  30  X  30  BAND  MATRIX. 
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TABLE  4.1:  REFERENCE  DATA  FOR  THE  WELLBORE  THERMAL  MODEL. 


Dr i 1 1  Pipe  ID  ( in. )  5 . 965 
Drill  Pipe  OD  (in.)  6.625 
Dri 1 1  Collar  ID  (in. )  2.750 
Dri 1 1  Collar  OD  (in. )  6.750 
First  Casing  ID  (in.)  19.124 
First  Casing  OD  (in.)  20.000 
Second  Casing  ID  (in.)  12.515 
Second  Casing  OD  (in.)  13.375 
Third  Casing  ID  (in.)  8.835 
Third  Casing  OD  (in.)  9.625 
Total  well  diameter  (in.)  26.000 
Open  hole  diameter  (in.)  8.375 

Total  Depth  (ft)  15000 

Depth  to  top  of  Drill  Collars  (ft)  13000 
Depth  of  First  Casing  Shoe  (ft)  2000 
Depth  of  Second  Casing  Shoe  (ft)  5000 
Depth  of  Third  Casing  Shoe  (ft)  10000 
Depth  to  Cement  -  First  Csg.  Annulus  (ft)  0 
Depth  to  Cement  -  Second  Csg.  Annulus  (ft)  1000 
Depth  to  Cement  -  Third  Csg.  Annulus  (ft)  4500 


DENSITY 
lb/cu  ft 

HEAT  CAP. 
Btu/lb-F 

THERMAL  C 
Btu/hr-f t 

Drill  Pipe 

487.00 

0.  1 

25.0 

Drill  Coll ars 

556.00 

0.  1 

25.0 

Casing 

503.00 

0.  1 

25.0 

Drilling  Fluid 

74.81 

0.4 

1  .0 

Cement 

196.00 

0.5 

0.4 

Ear  th 

165.00 

0.2 

1.3 

Flow  rate  (cu  ft/hr) 

=  1685 

Surface  temperature 

(deg.F)  = 

59. 

5 

Geothermal  Gradient 

(deg.F/ft)  = 

0. 

0127 

, 
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cuttings  to  the  surface  and  controlling  subsurface 
pressures.  It  is  here  that  a  fundamental  question  arises  - 
how  to  best  model  the  drilling  fluid  rheology?  Drilling 
fluids  have  progressed  from  being  little  more  than  clay 
suspensions  to  highly  complex  substances  both  Theologically 
and  chemically.  This  is  further  compounded  by  the  fact  that 
the  rheology  and  chemistry  can  vary  enormously,  even  during 
the  course  of  drilling  a  single  well,  depending  on 
prevailing  circumstances.  Obviously  some  form  of  consistency 
must  be  utilized,  at  the  cost  of  accuracy,  as  it  is 
completely  unfeasible  to  design  a  drilling  fluid  programme 
incorporating  all  the  vagaries  of  the  rheology  and  chemistry 
of  the  mud. 

The  only  valid  generalization  about  drilling  fluids  is 
that  they  are  non-Newtonian.  Even  so,  several  early  workers 
on  the  subject  e.g.  van  Olphen  (1950)  and  Cardwell  (1953) 
assumed  them  to  be  Newtonian,  probably  due  to  the  fact  that 
very  little  was  Known  of  non-Newtonian  fluid  behaviour. 
Unlike  Newtonian  fluids, which  have  the  constitutive  equation 


no  single  constitutive  equation  exists  to  describe  exactly 
the  relationship  between  the  shear  stress  and  shear  rate  of 
all  non-Newtonian  materials  over  all  ranges  of  shear  rates. 
Even  if  such  an  equation  could  be  developed  its  intricacy 
would  defy  engineering  application.  Slawomirski  (1975)  has 
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derived  a  constitutive  equation  for  time  independent 
drilling  fluids  which  illustrates  this.  Although  three  major 
categories  of  non-Newtonian  systems  are  recognized,  namely, 
time  independent,  time  dependent  and  viscoelastic,  only  the 
time  independent  system  has  received  any  degree  of  study. 
Fortunately  the  large  majority  of  industrial  non-Newtonian 
fluids,  including  drilling  fluids,  fall  into  this  category. 
The  time  independent  fluids  can  be  further  subdivided  into 
Bingham  plastic,  pseudoplastic  and  dilatant  fluids.  Numerous 
simplified  empirical  "models"  have  been  developed  to  relate 
the  shear  stress  to  the  shear  rate  for  these  fluids, 
especially  the  pseudoplastics  which  constitute  the  largest 
and  probably  the  most  important  class  of  non-Newtonian 
fluids.  Skelland  (1967)  has  summarized  the  most  important  of 
these  equations  and  the  diversity  of  drilling  fluids  is  such 
that  a  particular  type  of  drilling  fluid  probably  exists 
that  would  be  described  by  each  of  these  models.  Slawomirski 
(1975)  contends  that  the  majority  of  drilling  fluids  are 
thixotropic,  but  the  equations  to  describe  such  behaviour 
are  so  complicated  as  to  be  inapplicable  to  engineering 
problems.  Hence  it  is  generally  accepted  that  drilling 
fluids  are  typified  by  either  the  Bingham  plastic 
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or  Ostwald-deWaele  power  law  models. 
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(4.15) 


This  simplified  linearized  version  of  the  former  was  very 
much  in  favour  until  the  past  decade  when  the  advent  of 
complicated  polymer  muds  and  the  expanding  use  of  oil  base 
muds  emphasized  the  need  for  the  power  law  model.  The  power 
law  model  is  easily  applied  and  hence  the  large  majority  of 
the  research  on  non-Newtonian  flow  utilizes  this  model  as 
best  typifying  pseudoplastics.  Bingham  plastic  fluids,  on 
the  other  hand,  are  found  only  rarely,  although  high  solids 
drilling  fluids  are  well  described  by  this  model. 

When  the  functions  of  a  drilling  fluid  are  considered 
it  is  obvious  that  a  pseudoplastic  would  be  the  most 
appropriate  type  of  fluid.  It  is  shear  thinning  so  that  at 
the  high  shearing  rates  present  at  the  bit,  the  pressure 
drop  is  minimized,  whereas  at  the  low  shear  rates  in  the 
annulus  the  viscosity  is  increased,  thus  enabling  the  large 
volume  of  cuttings  to  be  efficiently  removed. 

However,  using  the  power  law  model  is  more  a  matter  of 
convenience  than  of  theoretical  validity,  as  it  has  certain 
disadvantages.  Drilling  fluids  typically  possess  a  yield 
value  which  cannot  be  accounted  for  by  this  model. 
Furthermore  the  power  law  model  predicts  infinite 
viscosities  and  zero  viscosities  in  the  limits  of  very  low 
and  very  high  shear  rates,  respectively.  Real  fluids, 
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however,  exhibit  a  finite  and  constant  viscosity  at  zero 
shear  rate.  The  use  of  this  model  also  requires  that  the  two 
defining  parameters  K  and  n'  remain  constant  over  the  entire 
range  of  shear  stress.  Fortunately  these  limitations  appear 
to  be  unimportant  for  drilling  applications,  if  the  drilling 
fluid  properties  are  assumed  independent  of  temperature. 

Nonetheless  some  have  considered  these  constraints  as 
justification  for  the  use  of  the  Bingham  plastic  model  and 
it  is  still  used  today  despite  its  limitations  being  much 
more  significant  than  those  of  the  power  law  model.  The 
Bingham  model  accounts  for  the  yield  values  typical  of  most 
drilling  fluids,  but  it  assumes  a  linear  relationship 
between  shear  stress  and  shear  rate  after  an  initial  yield, 
something  that  is  not  true  of  drilling  fluids.  Another 
negative  feature  is  that  no  explicit  relationship  can  be 
derived  between  the  shear  stress  and  the  volumetric  flow 
rate.  Furthermore  the  form  used  is  a  simplified,  linearized 
version  first  proposed  by  Caldwell  and  Babbit  (1941)  and  the 
simplification  has  been  shown  to  be  erroneous  under  certain 
circumstances  by  Hanks  and  Pratt  (1967).  However,  the 
simplified  model  is  still  used,  particularly  in  drilling 
hydraulics,  without  any  regard  to  the  limitations. 

Recently  attempts  have  been  made  to  resolve  the 
controversy  by  proposing  models  that  combine  the  obvious 
shear  thinning  features  of  drilling  fluids  with  a  yield 
stress.  The  Herschel -Bu lkley  model 
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has  been  considered  but  its  application  at  the  wellsite  is 
difficult  due  to  the  intricacy  of  the  equations  describing 
the  flow  behaviour  and  the  fact  that  no  explicit 
relationship  exists  between  the  shear  stress  and  the 
volumetric  flow  rate  for  this  model.  Zamora  and  Bleier 
(1977)  have  advocated  the  use  of  a  simplified  version  of 
this  model  and  Robertson  and  Stiff  (1976)  derived  an 
original  model  for  a  yield-pseudoplastic 
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that  has  certain  advantages  over  the  Herschel -BulKley  model. 

Whilst  recognizing  the  attractiveness  of  the 
yield-pseudoplastic  category  for  drilling  fluids,  little 
work  has  been  done  on  the  equations  necessary  to  describe 
the  behaviour  of  such  a  fluid  in  the  wellbore.  For  this  work 
it  is  proposed  to  use  the  power  law  model  as  being  the  best 
overall  compromise,  and  primarily  because  the  vast  majority 
of  the  relevant  theory  in  the  literature  is  based  upon  this 
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4.5  Non-Newtonian  Convective  Heat  Transfer  Coefficients 

The  non-Newtonian  behaviour  of  drilling  fluids  is 
something  that  all  previous  publications  on  the  subject  of 
wellbore  temperature  distribution  have  ignored.  No  mention 
is  made  of  any  investigations  to  ascertain  the  effect  of  the 
pseudoplasticity  of  a  drilling  fluid  on  the  convective  heat 
transfer  coefficients. 

All  previous  work  on  this  subject  has  used  the 
conventional  Sieder-Tate  correlation  to  estimate  the 
convective  heat  transfer  coefficients,  including  Raymond 
(1969)  and  Keller  et  al.  (1973).  Sump  and  Williams  (1973), 
whilst  recognizing  that  the  use  of  the  Sieder-Tate 
correlations  resulted  in  anomalously  low  temperatures , 
modified  this  correlation  by  regressing  on  temperature  data 
from  six  Gulf  Coast  wells.  Their  empirical  relationship 
resulted  in  modified  heat  transfer  coefficients  which  they 
claim  give  improved  temperature  profiles.  However,  the  fact 
that  data  from  only  six  wells  were  used  makes  their 
relationship  highly  suspect  for  general  use. 

Inherent  in  the  use  of  the  Sieder-Tate  correlation  is 
the  assumption  that  flow  is  turbulent.  Whilst  this  is 
generally  the  case  within  the  drill  string,  it  is  seldom 
true  in  the  annulus.  Thus,  the  use  of  the  Sieder-Tate 
correlation,  even  with  the  assumption  of  non-Newtonian  flow, 
is  invalid  for  determining  annulus  convective  heat  transfer 
coefficients . 

Given  the  shortcomings  of  the  previous  work  on  wellbore 
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temperature  distributions  with  regard  to  the  convective  heat 
transfer  coefficients  used,  an  investigation  was  carried  out 
to  compare  the  heat  transfer  coefficients  obtained  from 
non-Newtonian  correlations  with  the  Sieder-Tate  values,  and 
to  determine  how  significant  is  the  difference  between  the 
resulting  coefficients. 

An  extensive  review  of  the  literature  on  the  subject  of 
non-Newtonian  convective  heat  transfer  coefficients  revealed 
a  remarkable  dearth  of  research  in  this  area.  Such 
relationships  that  do  exist  are  either  purely  empirical  and 
relatively  simple,  or  analytical  and  invariably  very 
complicated.  The  latter  usually  allow  for  the  noni sotherma 1 
nature  of  such  parameters  as  heat  capacity,  but  their 
intricacy  makes  their  use  impractical  for  the  application 
envisaged  in  this  study.  However,  the  empirical  nature  of 
the  simple  relationships,  which  would  be  ideally  suited  to 
this  thermal  model,  means  that  each  relationship  gives  a 
different  result,  and  no  independent  criterion  exists  to 
determine  which  is  the  most  accurate. 

The  convective  heat  transfer  coefficients  obtained  from 
the  correlations  that  would  be  of  use  in  this  study  are 
summarized  in  Tables  4-2,  4-3  and  4-4  for  turbulent  pipe 
flow,  laminar  annular  flow  and  laminar  pipe  flow, 
respectively.  These  are  all  empirical  or  semi -empi r i cal 
relationships,  and  none  of  the  analytical  relationships  have 
been  i ncluded . 

It  is  evident  from  Table  4-2  that  the  Lakshmi narayanan 
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TABLE  4.2:  COMPARISON  OF  NON-NEWTONIAN  CONVECTIVE  HEAT 
TRANSFER  COEFFICIENTS  FOR  TURBULENT  PIPE  FLOW  OF  POWER  LAW 
FLUIDS. 


(1)  Lakshmi narayanan  et  al .  (1976) 

(2)  Clapp  (1961) 

(3)  Non-Newtonian  Sieder-Tate 

(4)  Metzner  and  Friend  (1959) 


h  ( Btu/hr -  ft 2  - °  F ) 


ri 

K 

Re 

Pr 

(1) 

(2) 

(3) 

(4) 

0.55 

0.1586 

2649 

48.7 

101.2 

81 .0 

108.5 

17.0 

0.60 

0. 1321 

2675 

48.3 

101.6 

84.1 

109.0 

22.1 

0.65 

0.1100 

2703 

47.8 

101  .9 

87.7 

109.6 

27.3 

0.70 

0.0916 

2733 

47.3 

102.3 

91.9 

110.1 

32.3 

0.75 

0.0763 

2764 

46.7 

102.7 

96.5 

110.7 

37.  1 

0.80 

0.0635 

2796 

46.2 

103.  1 

101.5 

111.3 

41.7 

0.85 

0.0529 

2829 

45.7 

103.5 

106.9 

111.9 

46.0 

0.90 

0.0441 

2863 

45.1 

103.9 

112.6 

112.5 

50.0 

0.95 

0.0367 

2898 

44.6 

104.4 

118.7 

113.2 

53.7 

■ 
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TABLE  4.3:  COMPARISON  OF  NON-NEWTONIAN  CONVECTIVE  HEAT 
TRANSFER  COEFFICIENTS  FOR  LAMINAR  ANNULAR  FLOW  OF  POWER  LAW 
FLUIDS  (INNER  AND  OUTER  WALLS  RESPECTIVELY). 

( 1 )  Skel 1  and  ( 1 967 ) 

(2)  Tanaka  and  Mitsuishi  (1975) 


n# 

NUSSELT 

(1) 

NO. 

(2) 

h  (Btu/hr- 
(1) 

f  t 2  -  ’  F  ) 

(2) 

0.55 

4.04 

3.53 

7.32 

5.79 

6.39 

5.05 

0.60 

4.03 

3.55 

7.29 

5.77 

6.43 

5.09 

0.65 

4.01 

3.58 

7.26 

5.75 

6.48 

5.13 

0.70 

4.00 

3.60 

7.24 

5.72 

6.52 

5.  16 

0.75 

3.98 

3.63 

7.21 

5.71 

6.57 

5.20 

0.80 

3.97 

3.65 

7.  19 

5.69 

6.61 

5.23 

0.85 

3.96 

3.68 

7.  17 

5.67 

6.66 

5.27 

0.90 

3.95 

3.70 

7.  16 

5.66 

6.70 

5.30 

0.95 

3.94 

3.72 

7.  14 

5.65 

6.74 

5.33 

Criterion  for  Nusselt  No.  (2)  :  >  9.1 


Annular  Length  =  1000  ft 
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TABLE  4.4:  COMPARISON  OF  NON-NEWTONIAN  NUSSELT  NUMBERS  FOR 
LAMINAR  PIPE  FLOW  OF  POWER  LAW  FLUIDS. 

( 1 )  Metzner  (1965) 

(2)  Charm  and  Merrill  (1959) 

(3)  Mizushina  and  Kuriwaki  (1968) 

(4)  Bird  (1959) 

( 5 )  Gr igul 1  (1 956 ) 

(6)  Skel 1  and  (1967) 


NUSSELT  NOS. 


n' 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

0.55 

6.88 

8.16 

72.51 

4.68 

4. 16 

4.28 

0.60 

6.81 

7.99 

64.71 

4.63 

4.18 

4.25 

0.65 

6.74 

7.86 

58.50 

4.58 

4.21 

4.23 

0.70 

6.69 

7.75 

53.47 

4.54 

4.23 

4.20 

0.75 

6.64 

7.66 

49.32 

4.50 

4.25 

4.  19 

0.80 

6.60 

7.59 

45.85 

4.47 

4.28 

4.17 

0.85 

6.56 

7.53 

42.92 

4.44 

4.30 

4.  15 

0.90 

6.52 

7.48 

40.41 

4.41 

4.32 

4.  14 

0.95 

6.49 

7.43 

38.24 

4.39 

4.34 

4.13 

Pipe  Length  =  1000  ft 
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et  al.  (1976)  and  the  non-Newtonian  Sieder-Tate  values 
compare  favourably  with  the  Newtonian  Sieder-Tate  value  of 
113  Btu/hr-f t 2- ° F  if  the  value  of  n'  is  assumed  to  be 
between  0.7  and  0.8,  a  typical  value.  The  correlation  of 
Clapp  (1961)  has  been  criticized  by  Metzner  (1965)  for 
experimental  errors,  which  probably  explains  the 
discrepancies  in  Table  4-2.  The  relationship  of  Metzner  and 
Friend  (1959)  is  invalid  in  this  particular  study  because 
the  necessary  criterion  (Equation  4.18)  is  not  satisfied. 


(4.18) 


(4.19) 


Table  4-3  indicates  a  favourable  comparison  between  the 
non-Newtonian  Nusselt  numbers  and  the  Newtonian  value  of 
4.12  for  laminar  flow  quoted  by  Keller  et  al.  (1973).  This 
is  also  the  case  in  Table  4-4  where  the  use  of  laminar  pipe 
flow  correlations  is  justified  by  neglecting  the  radial 
temperature  gradient.  Thus,  the  Nusselt  number  will  be  the 
same  for  both  the  inner  and  outer  annular  walls.  Several  of 
these  correlations,  namely  Tanaka  and  Mitsuishi  (1975), 
Metzner  (1965)  and  Charm  and  Merrill  (1959),  involve  the 
Graetz  number  which  means  that  the  Nusselt  number  becomes 
length  dependent.  Since  the  typical  dimensions  of  a  wellbore 
ensure  that  the  flow  is  fully  developed,  such  correlations 
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are  unsuitable.  The  di screpencies  in  the  values  determined 
from  the  correlation  of  Mizushina  and  Kuriwaki  (1968)  are 
propably  due  to  the  fact  that  they  based  their  experiments 
on  CIVIC  solutions,  which  are  viscoelastic,  not  pseudoplastic. 

As  indicated  in  chapter  5,  variations  in  the  convective 
heat  transfer  coefficients  have  no  significant  effect  on  the 
annular  temperature  profile.  Thus,  the  use  of  Newtonian  heat 
transfer  coefficients  is  valid.  However,  the  relationship  of 
Lakshmi narayanan  et  al.  (1976)  was  used  in  this  study  for 
theoretical  completeness  and  because  it  is  the  most  recent 
correlation.  This  relationship  is  given  in  Equation  4.20. 

St  =  0.0107Re'0-33  Pr‘0-67  (4.20) 

The  value  of  4.12  was  used  for  annular  Nusselt  numbers. 


4.6  Calculation  of  the  Energy  Source  Terms 

The  heat  generation  within  the  system  is  the  result  of 
several  different  energy  sources  namely: 

1.  Rotational  energy  due  to  the  work  required  to  rotate  the 
drill  string. 

2.  The  work  done  by  the  drill  bit. 

3.  Viscous  energy  due  to  the  friction  losses  inside  the 
drill  pipe,  drill  bit  and  annulus. 

The  calculation  of  the  thermal  energy  due  to  the  first  two 
sources  is  achieved  by  assuming  that  all  the  mechanical 
energy  input  by  the  rotary  drive  is  converted  to  thermal 
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energy  along  the  length  of  the  drill  pipe.  Keller  et  al. 
(1973)  divided  this  energy  60%  to  40%  between  the  drill  pipe 
and  the  drill  bit,  respectively. 

Thus,  these  terms  will  be  dependent  upon  each 
individual  drilling  rig,  and  would  be  calculated  by  a  simple 
relationship  between  the  rotating  speed  of  the  string,  and 
the  rotary  horsepower  required  to  maintain  this  speed.  The 
horse  power  would  be  converted  directly  to  thermal  energy. 

The  determination  of  the  energy  terms,  resulting  from 
the  frictional  presure  losses  within  the  system,  is  more 
difficult.  To  calculate  the  energy  input  due  to  the 
frictional  pressure  losses,  it  is  necessary  to  first 
determine  these  losses  in  the  drill  string,  drill  bit  and 
annulus,  which  in  turn,  requires  relationships  between  the 
pressure  drop  and  the  Known  flow  parameters  and  geometry. 

It  has  already  been  shown  that  the  rheology  of  drilling 
fluids  is  adequately  represented  by  the  Ostwald  deWaele 
power  law  model.  Therefore,  the  appropriate  relations  for 
determining  the  frictional  pressure  losses  for  power  law 
fluids  are  used  in  the  model.  Unfortunately,  there  are  a 
number  of  such  relations  in  the  literature,  and  it  is 
necessary  to  select  the  ones  most  applicable  to  this 
particular  model . 

An  assumption  already  made  is  that  flow  is  turbulent  in 
the  drill  pipe  and  laminar  in  the  annulus.  A  method  is 
included  in  the  model  to  verify  that  this  assumption  is 
valid  for  the  flow  conditions  pertaining  at  any  particular 
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time.  The  traditional  criterion  for  establishing  this  is  to 
compare  the  Reynolds  number  with  a  value  of  2100.  The 
Reynolds  number  used  is  that  derived  by  Metzner  and  Reed 
(1955)  for  time  independent,  non-Newtonian  fluids. 


(4.21) 


This  may  be  readily  modified  to  apply  to  pseudoplastic 
fluids  on  the  assumption  that  n'  and  K  remain  constant. 


(4.22) 


Re 


However,  the  selection  of  a  Reynolds  number  of  2100  as  the 
transition  value  has  been  criticized  by  numerous  authors. 
Dodge  and  Metzner  (1959)  have  indicated  that  the  value  is  a 
function  of  n'  and  of  temperature.  Ryan  and  Johnson  (1959) 
introduced  a  more  useful  definition  of  a  transition 
criterion,  which  was  subsequently  modified  by  Hanks  and 
Christiansen  (1962)  and  Hanks  (1963),  to  enable  the 
criterion  to  be  applied  to  power  law  pseudoplastics.  This 
criterion  is  a  function  of  n'  only  and  is  thus  independent 
of  temperature. 
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(4.23) 


The  use  of  Equation  4.23  as  a  transition  criterion  in  this 
model  has  certain  advantages  over  the  more  conventional 
criterion  of  comparing  the  Reynolds  number  to  2100.  The 
value  can  be  determined  with  reasonable  accuracy  and  can  be 
assumed  to  remain  constant,  provided  the  drilling  fluid 
properties  remain  unchanged. 

Different  relationships  are  required  to  determine  the 
frictional  pressure  losses  in  the  drill  pipe,  drill  bit  and 
annulus,  due  to  their  different  flow  behaviour  and  geometry. 

4.6.1  Drill  String 

Central  to  any  calculation  of  pressure  drop  within  a 
pipe  is  the  determination  of  the  Fanning  friction  factor.  In 
the  case  of  a  drill  pipe,  this  is  complicated  by  the  fact 
that  the  interior  wall  is  rough,  to  a  varying  degree. 
Unfortunately,  virtually  no  experimental  or  theoretical  work 
has  been  done  on  the  subject  of  non-Newtonian  flow  in  rough 
pipes  and  either  Newtonian  concepts  have  been  used  or  the 
pipe  is  considered  smooth.  If  the  former  is  assumed,  an 
estimate  of  the  roughness  height  is  required. 

Whereas  a  number  of  equations  exist  in  the  literature 
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to  determine 
flow  of  power 
has  developed 
However,  its 
experimental 


the  Fanning  friction  factor  for  the  turbulent 
law  fluids  in  a  smooth  pipe,  Torrance  (1963) 
the  only  one  applicable  to  rough  pipes, 
origin  is  purely  theoretical  and  no 
work  has  been  done  to  verify  it.  This  equation 


/f 


=  ^  109  7  +  6-0-^ 


(4.24) 


is  of  similar  form  to  the  familiar  Nikuradse  and  von  Karman 
equations  for  Newtonian  fluids.  The  one  problem  with  this 
equation  is  that  it  requires  a  value  for  the  roughness 
height  which,  for  something  like  rusty  commercial  steel 
pipe,  would  vary  enormously  and  be  very  difficult  to 
evaluate.  Fig.  4-6  indicates  that  for  moderately 
pseudoplastic  fluids  (n  >  0.7)  the  variations  in  the 
friction  factor,  calculated  by  Torrance's  equation,  with  the 
roughness  height  is  small  and  well  within  the  accuracy  of 
the  other  flow  parameters.  Hence,  the  roughness  height 
itself  is  not  a  critical  parameter  and  an  average  value 
would  suffice.  Usually  a  roughness  height  of  0.0006  inches 
is  used  for  pipe  flow  calculations,  although  this  figure  is 
suspect  for  old  pipe.  However,  on  the  basis  of  Fig.  4-6.  for 
moderately  pseudoplastic  fluids,  this  figure  would  be 
satisfactory . 

If  the  laminar  sub-layer  is  of  sufficient  thickness 
then  the  assumption  of  smooth  walls  is  valid.  Dodge  and 
Metzner  (1959)  derived  an  equation  for  the  friction  factor 
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as  a  function  of  the  generalized  parameters  first  used  by 
Metzner  and  Reed  (1959),  which  is  directly  applicable  to 
power  law  fluids. 


0.4 


(4.25) 


1.2 


n 


Since  then  various  relationships  have  been  derived, 
including  one  by  Chen  (1961)  who  modified  the  original 
equation  of  Dodge  and  Metzner  so  as  to  be  explicit  in  f. 
Solutions  of  a  number  of  different  equations  for  the 
friction  factor  in  smooth  pipes  are  summarized  in  Table  4-5, 
with  the  theoretical  values  compared  to  the  experimental 
values  of  Bogue  (1960).  The  relationships  used  are  those  of 
Dodge  and  Metzner  (1959),  Chen  (1961),  Mohammed  (1975), 
Torrance  (1963)  and  Tomita  (1961),  all  of  which  are  readily 
solvable  for  f.  Hanks  and  Ricks  (1975)  have  recently 
criticized  Dodge  and  Metzner' s  equation  and  propose  what  is 
probably  the  best  means  of  determining  the  friction  factor 
in  smooth  pipe  in  view  of  the  fact  that  only  the  value  of  n' 
need  be  known,  thus  ensuring  that  the  value  is  temperature 
independent.  Unfortunately  the  method  of  solution  is  rather 
laborious  and  more  appropriate  for  design  purposes.  Whereas 
Table  4-5  is  far  from  being  statistically  conclusive,  it 
does  indicate  the  validity  of  the  Dodge-Metzner 
relationship. 
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TABLE  4.5:  COMPARISON  OF  FANNING  FRICTION  FACTOR  EQUATIONS 
FOR  TURBULENT  FLOW  IN  SMOOTH  PIPES  WITH  THE  EXPERIMENTAL 
VALUES  OF  BOGUE. 

( 1 )  Bogue  (  1 960 ) 

(2)  Dodge  and  Metzner  (1959) 

(3)  Chen  ( 1961 ) 

(4)  Mohammed  (1975) 

( 5 )  Torrance  (  1 963 ) 

( 6 )  Tomi ta  (1961) 

FANNING  FRICTION  FACTORS 

n'  Re  (1)  (2)  (3)  (4)  (5)  (6) 

.895  98830  0.00414  0.00414  0.00400  0.00372  0.00391  0.00448 

.700  12137  0.00547  0.00580  0.00587  0.00486  0.00540  0.00717 

.445  12137  0.00475  0.00477  0.00482  0.00371  0.00402  0.00767 

.530  17374  0.00449  0.00430  0.00439  0.00349  0.00364  0.00642 

.465  12238  0.00435  0.00435  0.00444  0.00340  0.00363  0.00695 

.827  31473  0.00507  0.00508  0.00509  0.00438  0.00481  0.00574 

.677  8873  0.00619  0.00620  0.00626  0.00503  0.00580  0.00778 
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The  determination  of  the  pressure  drop  is  readily 
achieved  from  the  Fanning  friction  factor  by  means  of 
Equation  (4.26) . 


AP 


2fpV2L 

Dgc 


(4.26) 


In  the  unlikely  event  that  flow  is  laminar,  the  standard 
pressure  drop  -  friction  factor  relationship  would  be 
util ized. 

f  ■  if  <4-27) 


Of  these  various  relationships,  Torrance's  equation  for  the 
friction  factor  of  rough  pipes  (Equation  4.24)  is  the  most 
suitable,  with  a  roughness  height  of  0.0006  inches,  in  view 
of  the  fact  that  it  is  simple,  allows  for  roughness  and  is 
independent  of  temperature. 

4.6.2  Drill  Bit 

Because  of  the  intricacy  of  the  flow  behaviour  through 
the  drill  bit,  the  methods  used  to  determine  the  fluid 
velocity  and  pressure  drop  through  the  bit  have  been 
simplified.  An  Equation  for  determining  the  pressure  drop 
through  the  drill  bit  can  be  derived  from  the  Bernoulli 
equation  for  flow  through  an  orifice.  Since  the  drill  pipe 
cross-sectional  area  is  very  much  larger  than  the  nozzle 
area,  the  orifice  coefficient  is  assumed  to  equal  the 


discharge  coefficient.  The  value  for  the  discharge 
coefficient  is  assumed  to  be  0.95  by  Schuh  (1964). 
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4.6.3  Annu 1  us 

The  annulus  region  is  complicated  by  a  number  of 
factors  not  relevant  to  flow  in  the  drill  pipe,  namely  the 
volume  of  rock  cuttings  present,  a  factor  that  is  often 
ignored  in  drilling  hydraulics  studies.  Additional 
complications  are  the  varying  annular  size  (Fig.  4-1)  which 
is  often  considered  uniform,  and  the  eccentricity  of  the 
rotating  pipe.  In  a  deep  well  on  land  with  all  the  casing 
landed  at  the  surface  i .e.  no  liners  present,  the  assumption 
of  uniform  annular  size  is  a  reasonably  valid  one.  It 
becomes  more  questionable  where  a  long  bottomhole  assembly 
is  used,  a  liner  is  present  or  in  offshore  wells,  where  the 
riser  and  choke  line  are  a  complicating  factor. 

A  number  of  methods  exist  in  the  literature  for 
incorporating  the  eccentricity  of  a  rotating  pipe  into 
annular  flow  equations.  However,  all  assume  that  the 
eccentricity  is  regular  and  the  degree  of  eccentricity  is 
known.  Neither  applies  to  this  particular  case  and  the 
irregularity  of  the  eccentricity  necessitates  that  it  be 
neglected  with  regard  to  the  derivation  of  the  relevant  flow 
equations . 
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As  mentioned  above,  the  flow  behaviour  in  the  annulus 
should  ideally  be  laminar  and  this  is  normally  the  case. 
However,  as  with  pipe  flow,  Hanks'  stability  criterion 
(Equation  4.23)  can  be  applied  to  confirm  that  this  is  the 
case  in  each  section  of  the  annulus.  The  fact  that  the  flow 
should  be  laminar  removes  the  significance  of  the  rough 
surfaces  (Bowen  1961),  so  the  pressure  drop  can  be 
determined  by  means  of  a  smooth  walled  annular  flow  equation 
of  which  a  number  exist. 

The  original  work  in  this  field  was  that  of  Fredrickson 
and  Bird  (1958)  which  is  still  considered  the  most  accurate 
method  of  determining  annulus  pressure  losses.  Unfortunately 
the  exact  solution  for  power  law  fluids  in  an  annulus  is 
intricate  unless  use  is  made  of  graphs.  For  this  reason 
several  attempts  have  been  made  to  simplify  the  procedure  or 
to  derive  the  solution  by  different  means.  One  common  way 
was  to  modify  pipe  flow  equations  by  means  of  the  hydraulic 
radius  concept  and  although  this  resulted  in  simple 
equations  the  validity  of  this  concept  has  been  subjected  to 
a  great  deal  of  criticism.  Another  more  accurate,  and  less 
controversial  method,  involved  approximating  the  annulus  by 
two  parallel  plates.  A  number  of  the  most  applicable  methods 
are  compared  with  the  Fredrickson  and  Bird  method  in  Table 
4-6  by  solving  each  using  example  data  from  Skelland  (1967). 
As  with  pipe  flow  equations  this  table  is  far  from 
conclusive  but  it  does  indicate  that  any  of  the  methods  are 


* 


‘ 


,  I ' 


44 


TABLE  4.6:  COMPARISON  OF  PRESSURE  DROPS  OBTAINED  USING 
VARIOUS  ANNULAR  FLOW  EQUATIONS  TO  CALCULATE  THE  PRESSURE 
GRADIENT  USING  AN  EXAMPLE  IN  SKELLAND  (1967),  PAGE  115. 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


Fredrickson  and  Bird  (1958) . 60.85 

Fredrickson  and  Bird  (  Parallel  Plate  )..61.50 
Fredrickson  and  Bird  (  Annular  Slit  )....54.90 

Savins  (  1958 ) . 61.20 

Mishra  and  Mishra  (1976) . 60.69 

Wallick  and  Savins  (1969) . 60.75 


Ib/f 1 2/f t 
lb/ft  2/ft 
lb/ft  2/f t 
lb/ft  2/ft 
lb/ft  2/ft 
lb/f 1 2/f t 
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applicable.  1  However,  those  of  Savins  (1958)  or  the 
parallel  plate  approximation  of  Fredrickson  and  Bird  (1958) 
are  the  easiest  to  apply,  the  remainder  involving  the  use  of 
iterative  or  numerical  methods  of  solution.  It  is  the 
parallel  plate  approximation  of  Fredrickson  and  Bird  (1958) 
that  is  used  in  this  study  in  view  of  the  fact  that  the 
geometry  of  the  system  being  modeled  is  such  that  the 
parallel  plate  approximation  is  valid  for  the  annulus. 
(Equation  4.29) 


AP 


2KL 


(r  -  r  )g 
v  o  a/ac 


2 (n 1 +1  )q 


irn'(r0  +  ra)(ro  -  r  )' 


n 


(4.29) 


Thus,  the  frictional  pressure  losses  may  be  determined 
for  the  drill  string,  drill  bit  and  annulus  using  equations 
4.26,  4.28  and  4.29,  respectively.  The  sum  of  these 
frictional  pressure  losses  should  equal  the  pump  circulating 
pressure  less  the  frictional  pressure  losses  in  the  surface 
connections.  From  the  pump  circulating  pressure,  flow  rate 
and  mechanical  efficiency,  the  input  horsepower  is 
calculated  by  Equation  4.30  and  divided  between  the  drill 
string,  drill  bit  and  annulus  by  the  ratios  of  their 
respective  frictional  pressure  losses. 


HP 


APqe 

1.98  x  106 


(4.30) 


1  The  annular  flow  equation  derived  by  Mishra  and  Mishra 
contains  a  typing  error  which  has  been  allowed  for  in  Table 
4-6.  Skelland' s  solution  is  wrong  and  has  been  corrected  by 
Mishra  and  Mishra. 
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These  values  are  converted  from  thermal  to  mechanical  energy 
by  Equation  4.31, 

E  =  2544.9  HP  (4.31 ) 

and  added  to  the  energy  values  already  determined  for  the 
work  done  at  the  drill  bit  and  the  work  required  to  rotate 
the  drill  string. 

The  resulting  values  are  the  energy  source  terms  used 
i n  the  model . 


. 
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5.  DISCUSSION  OF  RESULTS 


Once  the  computer  model  had  been  designed  and  made 
operational,  it  was  used  to  simulate  wellbore  temperature 
distributions  under  a  variety  of  conditions.  Two  versions  of 
the  computer  program  were  used,  namely; 

1.  A  "full  scale"  model  incorporating  a  full  well 
completion  of  up  to  three  casing  strings  and  using  the 
data  in  Table  4-1  to  give  a  reference  temperature 
distribution. 

2.  A  "simplified"  model  which  neglected  the  well  completion 
and  assumed  a  uniform  hydraulic  radius  with  neither 
casing  nor  cement  present.  The  relevant  data  in  Table 
4-1  were  also  used  by  this  model  to  produce  a  second 
reference  temperature  distribution. 

The  annulus  temperature  profiles  for  the  two  different 
versions  are  shown  in  Fig.  5-1. 


5.1  Comparison  with  Previous  Models 

It  is  not  possible  to  compare  the  annulus  temperature 
profiles  generated  by  this  model,  with  those  for  most 
previous  models  because  the  respective  authors  have  not 
included  sufficient  data.  Often  values  are  not  given  for 
parameters  which  are  shown  later  in  this  chapter  to  be  very 
significant,  e.g.  Raymond  (1969)  includes  all  data  except 
the  fluid  heat  capacity,  which  is  shown  to  be  the  most 
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significant  parameter  of  all  in  terms  of  its  influence  on 
the  temperature  profile. 

The  only  models  which  allow  direct  comparison  are  those 

• 

of  Holmes  and  Swift  (1970)  and  Keller  et  al.  (1973).  Since 
this  model  is  essentially  the  Keller  et  al.  model  but  with  a 
much  improved  solution  procedure,  the  annulus  temperature 
profiles  are  very  similar.  Fig.  5-1  shows  a  comparison 
between  the  temperature  profile  of  Holmes  and  Swift  and  the 
two  versions  of  this  model.  It  is  evident  that  the  steady 
state  model  of  Holmes  and  Swift  is  in  poor  agreement  with 
the  unsteady  state  model  developed  by  Keller  et  al.  and 
extended  in  this  work.  It  is  further  evident  that  the  added 
improvement,  of  considering  casing  strings  and  alterations 
in  the  hydraulic  radius,  has  only  a  marginal  effect  on  the 
annulus  temperature  profile.  The  annulus  temperature 
profiles  generated  by  the  two  versions  differ  by  a  maximum 
temperature  of  5CF,  and  the  bottomhole  temperatures  are 
almost  equal . 

The  annulus  temperature  profiles  for  the  two  versions 
of  this  model  in  Fig.  5-1  do  not  include  heat  generation 
within  the  system  in  view  of  the  fact  that  the  energy 
sources  were  ignored  by  Holmes  and  Swift.  If  heat  generation 
is  included  in  the  two  unsteady  state  temperature  profiles, 
as  Keller  et  al .  appear  to  have  done  in  their  comparison 
(Keller  et  al.  Fig.  2),  the  bottomhole  temperatures  are 
almost  identical  for  the  data  in  Table  4-1.  However,  in  this 
study  (see  Fig.  5-1)  it  was  found  that  the  Holmes  and  Swift 


. 
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annulus  temperature  profile  continuously  diverges  from  the 
unsteady  state  temperature  profile  up  the  annulus,  and  the 
outlet  temperatures  differ  by  about  20° F.  This  contradicts 
Keller  et  al.  Fig.  2  in  their  work  indicates  similar  outlet 
temperatures  for  the  unsteady  state  temperature  profile  and 
the  Holmes  and  Swift  temperature  profile.  Whether  energy 
sources  are  included  or  not,  agreement  between  the  steady 
state  Holmes  and  Swift  annulus  temperature  profile  and  the 
unsteady  state  annulus  temperature  profiles  of  this  model  is 
very  poor.  This  is  not  altogether  surprising  considering 
that  Holmes  and  Swift  used  different  and  simpler  equations 
to  describe  the  wellbore  heat  flow,  and  consequently  were 
able  to  use  an  analytical  solution. 

It  should  be  noted  that  where  a  long  string  of  drill 
collars  is  present,  this  must  be  allowed  for  in  calculating 
the  frictional  pressure  losses  for  the  drill  string  and 
annulus  in  the  "simplified"  model.  Large  errors  may  occur  in 
the  respective  energy  source  terms  if  only  the  drill  pipe 
were  considered. 

The  effects  of  considering  the  energy  sources  in  the 
system  are  shown  in  Fig.  5-2  and  agree  closely  with  a 
similar  graph  in  the  paper  of  Keller  et  al.  Fig.  5-2 
confirms  the  necessity  of  including  energy  sources  in  any 
model  that  attempts  to  estimate  such  temperature 
distributions,  in  view  of  the  fact  that  there  is  an  error  of 
about  30 °  F  in  the  bottomhole  temperature  when  heat 
generation  is  ignored. 
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5.2  Parametric  Sensitivity  Analysis 

The  major  advantage  in  carrying  out  such  an  analysis  is 
to  ascertain  to  what  extent  variations  in  a  particular 
parameter  will  affect  the  temperature  distribution  in  the 
wellbore.  This  is  of  special  importance  where  the  particular 
parameter  is  one  that  cannot  be  evaluated  directly,  and 
instead  must  be  given  an  assumed  value.  In  such  a  case,  the 
degree  to  which  variations  in  this  parameter  would  affect 
the  temperature  indicate  how  important  is  an  accurate  value 
for  the  parameter.  It  also  provides  a  qualitative  estimate 
of  the  error  inherent  in  the  model  caused  by  using  assumed 
values.  Such  parameters  would  be  those  relating  to  the 
earth,  such  as  the  geothermal  gradient. 

Some  parameters  which  have  a  strong  influence  on  the 
temperature  distribution  can  be  directly  and  accurately 
evaluated,  e.g.  circulating  time,  depth,  geometry  and 
drilling  fluid  characteristics.  Such  an  analysis  is 
therefore  less  important  with  regard  to  these  parameters, 
but  nonetheless  pertinent  in  view  of  the  fact  that  it  can  be 
used  to  judge  the  validity  of  assuming  some  of  them  to 
remain  constant  e.g.  the  drilling  fluid  properties. 

Any  model  of  such  an  imprecise  and  complicated  an 
operation  as  drilling  a  well  must,  by  necessity,  be  a 
simplification.  The  fewer  the  simplifications  and 
assumptions,  the  more  accurate  the  model.  Unfortunately, 
there  is  no  physical  means  of  measuring  the  dynamic 
temperature  profile  in  a  well,  and  hence  there  is  no 
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independent  way  of  assessing  the  accuracy  of  the  model.  For 
this  reason,  any  variation  in  a  parameter  that  caused  a 
maximum  deviation  in  the  annulus  temperature  profile  of  less 
than  5°F  was  considered  insignificant.  It  should  be  noted 
that  interaction  between  parameters  was  ignored  and  that  if 
several  of  the  parameters  were  varied  simultaneously,  it  is 
possible  that,  whereas  individually  such  variation  had  an 
insignificant  effect,  the  sum  effect  could  be  significant. 

With  regard  to  the  actual  design  of  the  computer 
program,  the  following  observations  were  made. 

1.  The  use  of  double  precision  arithmetic  instead  of  single 
precision,  had  no  effect  on  the  resulting  annulus 
temperature  profiles,  thus  permitting  a  100%  reduction 
in  the  computer  memory  storage  requirements. 

2.  The  number  of  time  steps  used  had  no  effect  on  the 
solution,  which  is  to  be  expected  from  a  direct, 
implicit  solution  method. 

3.  The  size  and  number  of  the  vertical  elements  in  the 
matrix  had  an  insignificant  effect  on  the  annulus 
temperature  profile  for  a  15,000  feet  well  as  long  as 
the  elements  were  less  than  500  feet  in  size.  It  was 
observed  that,  in  general,  the  vertical  element  size  had 
to  be  less  than  3%  of  the  total  well  depth  to  ensure 
that  the  annulus  temperature  profile  remained 
independent  of  the  vertical  element  size. 

4.  The  optimum  size  and  number  of  the  radial  elements  in 
the  matrix  was  more  difficult  to  assess.  If  the  radius 
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of  the  external  radial  boundary  is  too  small  then  there 
will  be  heat  buildup  along  this  boundary  because  of  the 
external  radial  boundary  condition,  which  prohibits  the 
flow  of  heat  across  the  boundary  into  the  surrounding 
formation.  Long  circulation  times  or  excessive  heat 
generation  will  cause  the  radial  temperature  gradient  to 
intersect  the  external  boundary  before  the  gradient  has 
achieved  a  constant  temperature.  Ten  radial  elements 
were  arbitrarily  used  by  Keller  et  al.  (1973),  of  which 
the  first  seven  were  spaced  in  such  a  way  as  to 
represent  the  various  individual  sections  of  the 
wellbore  (Fig.  4-2)  and  the  remaining  three  used  to 
represent  the  earth.  The  same  arrangement  was  used  in 
this  work  with  the  three  outer  radial  elements  being 
each  two  feet  wide,  so  as  to  satisfy  the  external  radial 
boundary  condition.  Varying  the  size  of  these  three 
outer  elements  within  one  order  of  magnitude  had  no 
significant  effect  on  the  annulus  temperature  profile 
provided  there  was  no  heat  buildup  along  the  boundary. 
However,  with  the  simplified  version  of  the  model  the 
radial  element  size  is  much  more  important  in  view  of 
the  fact  that  only  the  first  three  elements  can  be 
spaced  to  coincide  with  the  physical  dimensions  of  the 
pipe  and  annulus.  Reducing  the  number  of  radial  elements 
had  an  adverse  effect  on  the  annulus  temperature 
profile,  as  did  alterations  in  their  size.  If  the  seven 
outer  elements  are  all  of  the  same  size,  the  optimum 
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size  is  the  minimum  to  satisfy  the  external  radial 
boundary  condition  for  a  given  circulation  time. 
Unfortunately,  increasing  the  circulation  time  will 
cause  a  heat  buildup  along  the  boundary  and  consequently 
create  discrepancies  in  the  annulus  temperature  profile. 
This  situation  was  avoided  by  reducing  the  size  of  the 
fourth  to  seventh  elements  and  fixing  the  size  of  the 
three  outermost  elements  at  2  feet  each.  This  closely 
approximates  the  real  case  where  casing  strings  are 
present.  It  was  found  that  the  closest  agreement  between 
the  two  versions  was  achieved  when  the  fourth  to  seventh 
elements  were  each  0.5  feet  in  width  and  the  three  outer 
elements  were  each  two  feet  in  width.  Not  only  was  the 
annulus  temperature  profile  insensitive  to  changes  in 
the  element  size  for  this  case,  but  no  heat  buildup 
occurred  along  the  external  boundary,  irrespective  of 
the  circulation  time. 

5.  The  steady  state  model  of  Holmes  and  Swift  has  already 
been  shown  to  be  inadequate.  A  more  meaningful 
comparison  of  steady  state  and  unsteady  state  annulus 
temperature  profiles  is  achieved  in  Fig.  5-3  where  an 
assumed  steady  state  annulus  temperature  profile  was 
generated  by  solving  the  model  for  a  circulation  time  of 
1,000  hours.  Not  only  is  this  profile  completely 
different  from  that  of  Holmes  and  Swift,  but  it  further 
supports  the  contention  that  the  steady  state 
approximation  of  wellbore  heat  transfer  is  not 
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justified.  Even  after  circulating  for  100  hours,  the 
maximum  feasible  time  of  continuous  circulation  likely 
for  any  bit  run,  the  bottomhole  temperature  differed 
from  the  steady  state  one  by  almost  20* F. 

An  anomaly  appears  to  exist  between  Fig.  2  and  Fig. 
3  in  the  work  of  Keller  et  al .  Whereas  Fig.  3  in  their 
work  agrees  closely  with  Fig.  5-3  in  this  study,  no 
agreement  exists  between  the  24  hour  temperature 
profiles  in  their  two  figures.  The  data  used  to  generate 
Fig.  2  and  Fig.  3  are  presumably  the  same  so  the  24  hour 
curve  in  Fig.  2  should  match  one  of  the  curves  in  Fig. 

3,  depending  on  whether  heat  generation  was  included  in 
the  generation  of  Fig.  2.  No  agreement  exists  and 
insufficient  detail  is  included  in  their  paper  to 
explain  the  anomaly.  The  unsteady  state  profiles  shown 
in  Fig.  2  of  Keller's  paper  could  not  be  reproduced  by 
the  model  developed  in  this  study,  although  the  curves 
in  Fig.  3  of  Keller's  paper  could  be  reproduced,  as 
shown  in  Fig.  5-3. 

The  following  parameters  were  found  to  have  a  significant 

influence  on  the  annulus  temperature  profile. 

1.  Drilling  fluid  heat  capacity 

Increasing  or  decreasing  this  parameter  by  50%,  which 
covers  the  range  of  values  quoted  in  the  literature, 
resulted  in  large  variations  in  the  annulus  temperature 
profile  as  shown  in  Fig.  5-4.  Since  a  small  change  in 
the  value  of  this  parameter  results  in  a  large 
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alteration  in  the  annulus  temperature  profile,  and  as 
the  variation  of  this  parameter  within  the  system  cannot 
be  accurately  assessed,  it  is  highly  significant. 

2.  Drilling  fluid  density 

Increasing  the  fluid  density  from  75  lb/ft3  to  100 
lb/ft3  decreased  the  bottomhole  temperature  by  22° F, 
whilst  increasing  it  to  120  lb/ft3  decreased  the 
bottomhole  temperature  by  38° F.  The  resulting  annulus 
temperature  profiles  are  shown  in  Fig.  5-5.  Although  the 
fluid  density  is  Known  with  reasonable  accuracy  upon 
entering  the  system,  its  variations  within  the  system 
are  not.  In  particular,  the  fluid  density  in  the  annulus 
will  increase  due  to  the  drill  cuttings.  The  magnitude 
of  this  increase  will  depend  on  the  drilling  rate,  the 
hydraulic  radius  and  the  circulation  rate  and  can  be 
determined  approximately  from  Equation  5.1. 


The  effect  of  this  increase  is  shown  in  Fig.  5-6  and 
obviously  should  not  be  ignored,  as  has  been  the  case 

with  all  previous  models. 

The  range  of  fluid  densities  considered  in  Figs. 
5.5  and  5.6  reflects  the  maximum  likely  fluctuation  in 
densities  due  to  gas  or  water  influx  to  the  annulus  or 
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DEN8ITY  =  76  LB8/CU.FT 
DENSITY  =  100  LB8/CU.FT 
DEN8ITY  =  120  LB8/CU.FT 
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drill  cuttings  buildup  in  the  annulus. 

3.  Fluid  flow  rate 

A  25%  decrease  in  the  flow  rate  resulted  in  a  25  #F 
increase  in  the  bottomhole  temperature.  Increasing  the 
flow  rate  by  25%  resulted  in  a  20' F  decrease  in  the 
bottomhole  temperature.  Although  the  flow  rate  has  a 
significant  effect  on  the  annulus  temperature  profile 
(Fig.  5-7  )  it  can  be  accurately  measured. 

4.  Geothermal  gradient 

Alterations  in  the  geothermal  gradient  have  a  marked 
effect  on  the  annulus  temperature  profile  as  shown  in 
Fig.  5-8.  Although  this  parameter  will  remain  constant 
for  a  particular  well,  it  cannot  be  known  with  any 
degree  of  accuracy.  Fig.  5-8  indicates  that  a  poor 
assumption  of  the  value  of  the  geothermal  gradient  will 
cause  errors  in  the  annulus  temperature  profile.  The 
range  of  values  in  Fig.  5-8  is  representative  of  the 
values  quoted  in  the  literature. 

5.  Earth  thermal  conductivity 

A  four-fold  variation  in  the  value  for  the  earth's 
thermal  conductivity  had  a  significant  effect  on  the 
annulus  temperature  profile  as  shown  in  Fig.  5-9. 
Reducing  it  by  400%  decreased  the  bottomhole  temperature 
by  20° F,  whilst  increasing  it  by  400%  had  only  a  slight 
effect  on  the  annulus  temperature  profile.  Although  a 
400%  alteration  in  this  parameter  within  a  wellbore  is 
highly  unlikely,  it  represents  the  complete  range  of 
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possible  thermal  conductivities.  It  is  a  parameter  that 
cannot  be  measured  directly  and  Fig.  5-9  indicates  that 
a  reasonable  estimate  of  the  earth's  thermal 
conductivity  is  required. 

6.  Hydraulic  radius 

Fig.  5-10  indicates  that  the  annulus  temperature  profile 
is  relatively  insensitive  to  this  parameter  with  a 
maximum  variation  in  the  annulus  temperature  profile  of 
about  1 0 *  F .  This  is  further  justification  for  the  use  of 
the  "simplified"  version  of  the  model  in  view  of  the 
fact  that  it  assumes  a  uniform  hydraulic  radius.  The 
variations  in  the  hydraulic  radius  in  Fig.  5-9  are  over 
the  length  of  the  well  and  the  assumption  of  a  constant 
value  would  involve  a  much  smaller  error  in  the  annulus 
temperature  profile. 

7.  Inlet  fluid  temperature 

This  can  be  accurately  measured  although  Fig.  5-11 
indicates  that  it  must  be  continuously  monitored  during 
the  operation  of  the  model  at  the  wellsite  to  ensure 
that  the  annulus  temperature  profiles  are  accurate. 
Reasonable  variations  in  the  earth  density  and  specific 
heat,  the  drilling  fluid  thermal  conductivity,  the 
convective  heat  transfer  coefficients  and  the  drill  pipe 
density  and  thermal  properties,  had  no  significant  effect  on 
the  annulus  temperature  profile. 
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FIG. 5-10:  EFFECT  OF  VARIATIONS  IN  THE  HYDRAULIC  RADIUS 

ON  THE  ANNULUS  TEMPERATURE  PROFILE. 
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5.3  Validity  of  the  Assumptions 

There  is  no  evidence  from  this  work  to  suggest  that  any 
of  the  major  assumptions  enumerated  in  Chapter  4  are 
invalid.  Given  the  physical  geometry  of  the  system  being 
modeled,  they  are  all  logical  assumptions  and  the  problem  is 
greatly  simplified  as  a  result. 

However,  the  final  assumption,  namely  that  the  fluid 
properties  are  independent  of  temperature,  is  open  to 
question.  It  is  evident  from  Figs.  5.4  and  5.5  that  the 
fluid  specific  heat  and  fluid  density  have  a  significant 
effect  on  the  temperature  distribution.  No  investigation  was 
carried  out  to  assess  the  extent  to  which  these  parameters 
are  affected  by  temperature  and  until  this  is  done,  the 
validity  of  this  assumption  must  remain  in  doubt. 

With  regard  to  the  assumptions  made  by  previous 
workers,  the  implied  assumption  of  Newtonian  flow  behaviour 
is  valid  from  both  the  point  of  view  of  convective  heat 
transfer  coefficients  and  the  calculation  of  the  energy 
source  terms.  Although  this  assumption  gives  increased 
values  for  the  frictional  pressure  in  the  drill  pipe  and  the 
annulus,  and  hence  energy  source  terms,  the  error  causes  an 
insignificant  alteration  in  the  annulus  temperature  profile. 

The  assumption  of  steady  state  heat  transfer  and 
ignoring  the  energy  sources  in  the  system,  have  already  been 
shown  to  be  unjustified. 
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5.4  Wellsite  Operation  of  the  Model 

Information  on  the  wellbore  temperature  distribution  is 
most  likely  to  be  required  after  the  surface  and  conductor 
casing  strings  have  been  set.  The  operation  of  this  model 
would  be  most  suited  to  periods  of  constant  drilling  through 
thick,  uniform  formations,  a  situation  most  likely  to  occur 
after  several  thousands  of  feet  have  been  drilled.  Because 
the  initial  temperatures  in  the  system  are  assumed  to  be  the 
geothermal  gradient,  this  will  only  be  valid  prior  to  the 
commencement  of  drilling  when  the  temperature  profile  will 
not  have  been  disturbed.  Thus,  for  the  operation  of  the 
model  to  begin  at  some  point  during  the  drilling  of  the 
well,  an  estimate  must  be  made  of  the  initial  temperatures 
prevailing  at  that  time. 

If  the  operation  of  this  computer  model  is  assumed  to 
commence  after  a  casing  string  has  been  set,  a  logical 
starting  point,  then  two  options  exist  for  determining  the 
initial  temperatures.  Either  the  static  bottomhole 
temperature  obtained  from  a  logging  tool  could  be  used  to 
give  a  "pseudo-geothermal  gradient",  or  if  sufficient  time 
had  elapsed  since  the  last  period  of  circulation,  the 
temperature  distribution  in  the  wellbore  would  have  reverted 
approximately  to  the  geothermal  gradient.  Fig.  5-12 
indicates  that  the  annulus  temperature  profile  reverts 
almost  to  the  geothermal  gradient  after  about  2  days  without 
circulation.  This  is  about  the  minimum  time  required  to  run 
a  casing  string. 
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For  the  static  case,  heat  transfer  is  by  radial  and 
vertical  conduction  only  and  convection  is  negligible.  Thus, 
the  third  assumption  in  Chapter  4  is  no  longer  valid  and 
Equations  4.1  and  4.3  should  be  modified  accordingly.  This 
was  not  done  and  instead  the  flow  rate,  heat  generation  and 
convective  heat  transfer  coefficients  were  reduced  to  zero. 
Hence,  Fig.  5-12  is  only  an  approximation. 

The  computer  program  should  be  initiated  when 
circulation  first  commences  after  a  casing  run  with  either 
the  geothermal  gradient  or  the  gradient  obtained  from  a 
logging  tool  bottomhole  temperature,  used  as  the  initial 
temperature  profile. 

Once  the  computer  program  has  commenced  generating 
temperature  data,  it  may  be  left  to  run  continuously, 
obtaining  the  necessary  data  from  electronic  sensors  after 
predetermined  time  intervals  have  elapsed.  Hence,  the 
salient  drilling  parameters  such  as  flow  rate,  drilling 
fluid  density,  and  drilling  fluid  inlet  temperature  can  be 
continuously  updated.  Should  a  means  be  available  to 
accurately  measure  the  drilling  fluid  heat  capacity  at  the 
wellsite,  this  parameter  should  also  be  updated.  In  fact  it 
is  an  essential  requirement  for  the  accurate  operation  of 
this  model.  If  drilling  ceases,  this  can  be  allowed  for  by 
considering  the  flow  rate,  heat  generation  and  the 
convective  heat  transfer  coefficients  as  being  equal  to 


zero. 


. 
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6.  SUMMARY  AND  CONCLUSIONS 


The  primary  objective  of  this  work  was  to  develop  a 
computer  model  that  estimates  the  temperature  distribution 
in  a  wellbore  both  accurately  and  quickly.  This  model  was 
then  used  to  examine  the  validity  of  the  assumptions  made  in 
this  and  previous  models,  and  to  investigate  the 
significance  of  the  parameters  in  terms  of  their  effect  on 
the  temperature  distribution.  The  results  of  this 
investigation  contradicted,  in  several  instances, 
conclusions  reached  by  previous  workers. 

The  major  conclusions  of  this  study  are  as  follows: 

1.  Using  a  band  algorithm,  direct  solution  technique  to 
solve  the  finite  difference  equations  describing  the 
transient  heat  flow  in  a  wellbore,  a  solution  time  of 
about  3  seconds  was  required  on  the  University  of 
Alberta's  Amdahl  470  computer  as  compared  to  170  seconds 
on  a  CDC  6400  quoted  by  Keller  et  al.  (1973).  The  Amdahl 
470  is  approximately  five  times  faster  than  the  CDC 
6400,  so  the  solution  method  used  in  this  work  is  an 
order  of  magnitude  faster  than  the  method  used  by  Keller 
et  al  . 

2.  The  following  parameters  have  a  significant  effect  on 
the  wellbore  temperature  distribution. 

a.  drilling  fluid  heat  capacity 

b.  drilling  fluid  density 

c.  drilling  fluid  flow  rate 
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d.  drilling  fluid  inlet  temperature 

e.  geothermal  gradient 

f.  thermal  conductivity  of  the  earth 

g.  time  of  circulation 

h.  depth 

3.  Energy  generation  within  the  system  must  be  allowed  for. 

4.  The  allowance  for  a  number  of  casing  strings  and  a 
string  of  drill  collars  is  unnecessary,  and  a  wellbore 
of  uniform  hydraulic  radius  is  adequate. 

5.  The  assumption  of  steady  state  wellbore  heat  transfer  is 
unwarranted . 

6.1  Suggestions  for  Further  Research 

Variations  in  the  drilling  fluid  density  and  heat 
capacity  have  a  significant  effect  on  the  wellbore 
temperature  distribution.  However,  one  of  the  assumptions 
made  in  this  study  is  that  the  fluid  properties  remain 
constant,  thus  simplifying  the  partial  differential 
equations  describing  the  wellbore  heat  transfer.  If  in  fact, 
these  two  parameters  are  not  temperature  invariant,  then  the 
temperature  distribution  will  be  affected,  and  errors  will 
occur . 

It  is  therefore  recommended  that  the 
temperature-dependent  behaviour  of  these  two  parameters  be 
investigated  to  ascertain  whether  the  increased  accuracy  of 
allowing  for  this  behaviour  justifies  the  increased 
intricacy  of  the  partial  differential  equations  and  their 
subsequent  solution. 


. 
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It  is  also  essential  that  a  means  be  available  at  the 
wellsite  to  accurately  measure  the  drilling  fluid  heat 
capacity.  This  would  constitute  part  of  the  standard 
procedure  for  measuring  the  fluid  properties.  An 
investigation  into  some  simple  method  of  obtaining  this 
information  would  be  of  value. 

There  are  no  means  available  to  measure  the  annulus 
temperature  profile  nor  the  dynamic  bottomhole  temperature, 
so  no  means  exist  for  determining  how  precisely  this,  or 
previous  models,  predict  actual  temperatures  in  the  field. 
The  only  temperatures  that  can  be  measured  in  the  field  are 
the  inlet  and  outlet  fluid  temperatures .  However,  in  view  of 
the  fact  that  the  model  will  estimate  the  annulus  surface 
temperatures,  some  idea  of  its  accuracy  could  be  obtained  by 
comparing  the  measured  annulus  outlet  temperatures  with 
those  obtained  from  the  model. 

The  model  also  allows  for  the  estimation  of 
temperatures  during  periods  of  no  flow,  so  the  bottomhole 
temperature  in  the  wellbore  after  a  period  of  circulation, 
which  may  be  obtained  from  a  logging  tool,  should  be 
estimated  by  the  model.  The  comparison  between  these  values 
will  provide  further  information  on  the  accuracy  of  the 
model . 
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APPENDIX  I 


DERIVATION  OF  THE  PARTIAL  DIFFERENTIAL 


EQUATIONS  DESCRIBING 


THE  TRANSIENT  HEAT  FLOW  IN  A  WELLBORE 
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The  differential  element  used  to  derive  the  partial 
differential  equations  describing  the  transient  heat  flow 
in  a  wellbore  is  depicted  below. 


The  energy  balance  on  the  incremental  volume  may  be 
expressed  as  follows  ^ 

Input  +  Generation  =  Output  +  Accumulation 


' 
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FLUID  IN  THE  DRILL  PIPE: 
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Dividing  through  by  az  and  taking  limits  as  az  0 
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But  from  Fourier's  Law,  q  =  -k 
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2r.h.(T.  -  T  ) 


3T. 


w  3z2  K  -  rp)  ' p 


FLUID  IN  THE  ANNULUS: 


w 


<1  -  y 


=  P...C 


W 


w  pw  at 


(A. 1.8) 


pqCpTa 


+  Q  AZ 
z+az  wa 


pqCpTa|z  +  2*raha(Ta  '  VAZ 


+  2,Troho(Ta  ■  Tf)AZ 


+  jit  Wro  -  ra)pCPTaAz] 


(A. 1.9) 


Rearranging 


QaAZ  =  -  pqCp(Ta 


z+Az  "  a 


z)  +  2,raha(Ta  -  TJaz 


+  2irro%^Ta  '  Tf)AZ  +  £  -  ra)pCPTaAz]  (AJJ0) 


Dividing  through  by  az 


Qa  =  -  pqCr 


(T. 


z+az  ~  ^a 


AZ 


+  2,rraha(Ta  ‘  Tw) 


+  2,rroho(Ta  *  V  +  it  W*?  ’  ra)pCpTa] 


(A. 1.11) 


‘ 
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Taking  limits  as  az  -*  0 


^cp  ST  +  2*raha<Tw  -  V  +  2*roho(Tf  -  Ta> 


=  pC^7r(rJ  -  rl) 


aT 


p  v  o  a'  at 


(A. 1.12) 


FORMATION: 

The  general  energy  equation  for  a  two-dimensional  system, 
neglecting  the  radial  temperature  gradient  and  viscous  dissipation, 
may  be  written  as  - 


aL 


aL 


pfCpf  3t  +  VZ  3Z 


r  ar 


(A. 1.13) 


V  =  0  in  the  formation 

•m 

Therefore 


ar 


1 

r 


aT, 


ar 


pfcpf  3Tf 
k-  3t 


(A. 1.14) 


. 
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DERIVATION  OF 


APPENDIX  II 
THE  FINITE  DIFFERENCE 


EQUATIONS 
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Writing  the  partial  differential  equations  derived  in 


Appendix  I  in  finite  difference  form  involves  representing  the  first 
order  derivatives  by  two-point  forward  and  backward  difference 
approximations,  and  the  second  order  derivatives  by  three-point  centred 
difference  approximations. 

The  four  partial  differential  equations  can  be  written  in 
implicit  finite  difference  form  for  each  vertical  element  j  of 
thickness  az.,  and  for  a  time  step  of  At  as  follows  - 

J 


FLUID  IN  THE  DRILL  PIPE: 


pqc 


p 


p 


(A. 2.1) 


DRILL  PIPE  WALL: 


At 


(A. 2. 2) 


' 

. 


87 


FLUID  IN  THE  ANNULUS 


qcr 


n+1  n+1 

3J+T"  T3 


H) 

ill  + 


Zj 


2irr  h 
a  a 


i+l 

2,j 


-  T 


.n+1 
3  s  j 


.  .  . Tn+1  n+1 

+  2irr  h  ( T/.  .  T o  . 
o  o  \  4,j  3,j 


'1*0 


rZMP 


(i?\  -  tIJ 
y  3  ,j _ 3  9  j 


At 


-  Q. 


(A. 2. 3) 


FORMATION: 


a2T 


f  _  1 


dZ‘ 


Azj 


n+1  n+1 

Ti,J+1  "  » > J, 


AZ 


j+0.5 


n+1  n+1 

i ,  j  ~  i  ,j-l 

AZj-0.5 


(A. 2.4) 


Le.  _l  ( r  !Lf 

r  ar  \  ar 


ri  i  Ar 


ri -0 . 5  3r  ^Ti-0.5^  '  ri+0.5  3r  ^Ti+0.5 


(A. 2. 5) 


3r  ^Ti -0 . 5 ^  Ar  ^Ti-1  '  Ti 


(A.2.6) 


3r  ^Ti+0.5^  ~  Ar  ^Ti  *  Ti+1^ 


(A. 2. 7) 


From  Fourier's  Law 


dT 


f  _ 


q  -  2urkf -gj—  -  2^1+0.5  kf 


(Vi  -  V 

(ri+l  *  ri> 


(A. 2. 8) 


88 


For  radial  heat  flow 


q  = 


2  MT1  -  T,) 
1  n  (rg/r-j ) 


2)  _  2  kf(T1+1  -  Tt) 


In 


(ft) 


Equating  Equations  A. 2. 8  and  A. 2. 9 


i+0.5 


ri+i  -  ri 


In 


ri+l 


ri  -0 . 5 


ri  ‘  ri-l 


In 


ft) 


Substituting  into  Equation  A. 2. 5 


i  . 

r  ar  V  ar  / 


(ft-i  -  V 


(A. 2. 9) 


(A. 2. 10) 


(A. 2. 11) 


(A. 2. 12) 


Therefore 


89 


At 


(A. 2. 13) 


where 

r.  =  2r  +  (2i  -  7)Ar 
1  o 

Ar  =  Width  of  each  of  the  external  radial  elements 


5  <  i  <10 
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APPENDIX  III 
SUBROUTINE  LEQT1B 


■ 

. 
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The  linear  system  of  equations  describing  the  transient 
heat  flow  in  the  wellbore  is  expressed  in  matrix  form  in 
Equation  A . 3 . 1 

[AMT}  =  {C}  (A. 3.1) 

where  A  is  the  coefficient  matrix  in  band  storage  form,  T  is 
the  unknown  temperature  vector  at  the  new  time  level  and  C 
is  the  vector  containing  the  right  hand  sides  of  the 
equations . 

Solution  of  this  equation  to  obtain  the  values  of  the  T 
vector  is  achieved  by  means  of  the  International  Mathematics 
and  Statistics  Library  (IMSL)  subroutine  LEQT1B.  This 
subroutine  solves  the  general  equation  A. 3.2  by  first 
factorizing  the  coefficient  matrix  into  lower  and  upper 
triangular  matrices. 

[ A ] . [ X ]  =  [B]  (A. 3.2) 

It  then  solves  for  each  element  of  the  unknown  matrix  X  by 
backward  substitution. 

In  the  specific  case  being  considered  here,  the 
matrices  X  and  B  are  replaced  by  vectors  T  and  C, 
respectively. 

The  sample  output  included  in  Appendix  IV  was  obtained 
from  the  data  of  Table  4-1.  Once  the  computer  program  had 
calculated  the  coefficients  in  matrix  A  and  the  right  hand 
terms  contained  in  vector  C,  LEQT1B  was  used  to  solve 
Equation  A. 3.2  for  T. 

The  parameter  list  of  subroutine  LEQT1B,  with  a  brief 
definition  of  each  of  the  parameters  and  the  corresponding 
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value  of  the 

parameter  for  this  particular 

example,  are 

given  below. 

Parameter 

Def i ni t ion 

Value 

A 

Coefficient  matrix 

N 

Order  of  matrix  A 

1560 

NLC 

Number  of  lower  codiagonals 

12 

NUC 

Number  of  upper  codiagonals 

10 

IA 

Number  of  rows  in  A 

1560 

B 

Right  hand  side  matrix 

M 

Number  of  columns  in  B 

1 

IB 

Number  of  rows  in  B 

1560 

I  JOB 

Operation  flag  (see  below) 

0 

XL 

Work  area  of  minimum  size  N 

x  (NLC  +  1 )  21000 

IER 

Error  parameter  in  output 

The  mesh-centred  grid  representing  the  15,000  feet 
wellbore  is  composed  of  156  vertical  elements  and  10  radial 
elements.  After  the  steps  outlined  in  Chapter  3  have  been 
carried  out,  this  is  transformed  into  a  1560  x  21  band 
storage  matrix.  The  addition  of  the  two  extra  columns  to 
allow  for  the  equalization  of  temperatures  at  the  drill  bit 
increases  the  band  storage  matrix  to  1560  x  23  elements. 
Thus,  the  coefficient  matrix  A,  in  band  storage  form, 
comprises  1560  rows,  12  columns  to  the  left  of  the  main 
diagonal  and  10  columns  to  the  right  of  the  main  diagonal. 


When  the  operation  flag  I  JOB  is  zero,  the  complete 
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solution  procedure  is  carried  out  with  the  solution  matrix 
replacing  the  right  hand  side  matrix  on  output.  In  this 
case,  the  solution  vector  T  replaces  the  right  hand  side 
vector  C  which  then  contains  the  new  temperatures. 

If  IJOB  =  1,  only  the  factorization  occurs  and  if  IJOB 
=  2,  only  the  backward  substitution  takes  place.  Hence,  IJOB 
=  2  presupposes  that  LEQT1B  has  already  been  called  with 
IJOB  =  1  or  0. 


